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The interplay of the tunneling transfer of charges and the emission and absorption of light can 
be investigated in a setup, where a voltage-biased Josephson junction is connected in series with a 
microwave cavity. We focus here on the emission processes of photons and analyze the underlying 
time-dependent statistics using the second-order correlation function g f2 \r) and the waiting-time 
distribution w(r). Both observables highlight the crossover from a coherent light source to a single¬ 
photon source. Due to the nonlinearity of the Josephson junction, tunneling Cooper pairs can 
create a great variety of non-classical states of light even at weak driving. Analytical results for the 
weak driving as well as the classical regime are complemented by a numerical treatment for the full 
nonlinear case. We also address the question of possible relations between < 7 ^ 2 '(t) and w(r) as well 
as the specific information which is provided by each of them. 
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I. INTRODUCTION 

New possibilities offered by sources of nonconventional 
light are of interest for a wide range of applications. 
These range from secure quantum communication or 
teleportation to measurement, starting as early as the 
laser with more recent proposals exploiting squeezing 
properties for enhanced sensitivity and to tackle noise 
limitations 1 . On the other hand, superconducting circuit 
quantum electrodynamics devices have demonstrated the 
ability to create on demand any arbitrary state of light 
in a microwave stripline cavity, as shown by mapping 
out the Wigner densities of Fock states and their super¬ 
position, and of cat states in a multi-shot creation and 
read-out cycle 2 . Here we will discuss a setup combin¬ 
ing these ideas in exploiting a superconducting device as 
a continuous source of microwave light with interesting 
quantum-statistical properties. 

Over the last years, a variety of different experimen¬ 
tal realizations 3-8 and theoretical investigations 8-18 have 
been dedicated to the continuous emission of photons into 
a cavity mode by the Cooper-pair (CP) current across 
a clc-biased Josephson junction (JJ). A microwave res¬ 
onator connected in series to the junction hereby works 
as a well-defined electromagnetic environment absorbing 
the energy provided by the bias to each tunneling CP by 
creating photons. The nonlinearity of the JJ translates 
into a complex photon creation process stretching from 
suppression or enhancement of higher cavity excitations 
to multi-photon resonances. The excited photons finally 
leak from the resonator so that their mean emission rate 
and spectrum can be measured. 6 Employing high-Q su¬ 
perconducting resonators, far-from-equilibrium states 3,7 , 
with high-photon occupations can be achieved; this can 
be used as a basis for analyzing the quantum-classical 
crossover and to investigate nonlinear effects, 12-14 , such 
as number squeezing. Highly attractive in this class of 
systems is the feature that measuring the emitted mi¬ 
crowave radiation can indirectly yield information about 
the CP current and its fluctuations and vice versa. 6 


Generally speaking, these hybrid circuits combine ob¬ 
servational tools, methods, and phenomena known from 
the fields of quantum optics and charge transfer physics 
(Josephson photonics). Ultimate goals are to design 
sources for the on-demand production of single photons 
and highly non-classical multi-photon states (entangled 
photons) in the microwave up to the low-frequency tera¬ 
hertz regime. 

In this paper, we concentrate on the bright (photonic) 
side of this setup and study the time-resolved statistics 
of photonic emission events by means of the second-order 
correlation function g^ (r) and the waiting-time distri¬ 
bution w(t). Both functions are well-known tools from 
quantum optics 19,20 to study correlations of photonic 
emission events but have recently also been addressed to 
electron statistics in transport systems 21-26 , where rela¬ 
tions to other common statistical properties of electronic 
transport, e.g., full counting statistics or the noise spec¬ 
trum, have been discussed 21,24 . 

To model the circuit of a voltage-biased JJ in series 
with a resonator, we employ effective Hamiltonians de¬ 
rived in Refs. 12 and 13 to describe the system close to 
the one- and two-photon resonances. In conjunction with 
a master equation of Lindblad form modeling the impact 
of the environment on the system, most importantly the 
photon leakage, we can then study photonic correlations 
via the two time-dependent observables g^ 2 \r) and w(r). 
While the full nonlinear quantum case requires generally 
numerical calculations, analytical expressions can be ob¬ 
tained in the weak-driving regime as well as in the classi¬ 
cal limiting case on the basis of a perturbation treatment. 

Crucially, the inherent nonlinearity of the JJ enters 
the effective Hamiltonians as a nonlinear driving term, 
not as a nonlinear potential. The impact of nonlineari¬ 
ties is thus not only determined by the strength of driving 
given by the Josephson energy Ej , but more importantly 
by tunneling matrix elements of the drive Hamiltonian 
between resonator states. These depend on a dimension¬ 
less parameter k characterizing the importance of charge 
quantization of the CP current by the ratio of the to- 
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tal charging energy to the energy of a resonator photon. 
In consequence, the JJ-resonator system can reduce to a 
number of simpler systems, from a driven harmonic os¬ 
cillator or a parametric amplifier, when the CP current 
becomes classical (re —> 0 ) to a two- or few-level system 
for special values of re. 

We will investigate, how the crossover between these 
special cases is reflected in the statistics of emitted light, 
changing from a coherent light source to a single-photon 
source, where </ 2 )(r = 0) = 0. Additionally, we examine 
possible relations between g^(r) and w(t) and use the 
specific physical information which is provided by each 
of them to reveal to what extent the system features a 
renewal property. This is found to be the case for the 
harmonic-oscillator limit (re —>■ 0 or weak driving) as well 
as the two-level system. 

The paper is organized as follows. In Sec. II, we present 
our theoretical model and introduce the observables for 
studying correlations, the second-order correlation func¬ 
tion g^(r) and waiting-time distribution w(t). Sec¬ 
tion III discusses various aspects of the time-dependent 
statistics of photon emission events by means of these two 
observables close to the one-photon resonance (Sec. Ill A) 
and later in (Sec. IIIC) for the two-photon resonance. 
Relations between g^ 2> (t) and w(t) and renewal proper¬ 
ties of the system are considered in Sec. IIIB. Finally, 
Sec. IV contains the conclusions and an outlook address¬ 
ing the open questions for further research. 


II. MODEL 

In the following, we briefly discuss our modeling of 
the JJ-resonator circuit, i.e., we introduce the model 
Hamiltonian and a quantum master equation in Lind- 
blad form describing the impact of the environment. On 
the basis of the master equation, we can calculate ar¬ 
bitrary time-independent correlation functions, for in¬ 
stance, the second-order correlation function, </ 2 )(r), 
and the waiting-time distribution, w(t), introduced in 
the last part of this section. 


A. Hamiltonian 

The Hamiltonian of a voltage-biased JJ-resonator cir¬ 
cuit 

H= k + {^) ^ 2 -£ jC os(r?)-2 eV eS N (1) 

is constituted from its two subunits: the resonator de¬ 
scribed as harmonic oscillator with mass m = ( h/2e) 2 C 
and frequency ujq = 1/y/LC and a JJ with Josephson 
energy Ej and phase g across the junction . 13 These two 
parts are coupled dynamically by means of an effective 
voltage V e ff = V — V-es, where V represents the external 
voltage and V[ e s the voltage drop at the resonator. Here, 
the number operator N counts the number of CPs which 


have transferred across the JJ. Both charge operator q 
and resonator phase </> as well as number operator N and 
junction phase g represent sets of conjugate variables, 
i.e., [q,<fi\ = — 2 ie and [At, 77 ] = — i, respectively. 

Applying a gauge transformation on this Hamilto¬ 
nian ( 1 ), mapping it to a frame rotating with the 
(ac-) Josephson frequency, wj = 2 eV/h, and performing 
a rotating-wave approximation 2, finally yields 12,13 


TT (p ) = /jA n - 


(~i) p £j* 
2 


(a t ) p + (-l) p a 1 ’ 


J p (V^ren) 
n p / 2 

( 2 ) 


close to the p-photon resonance, A = wo — wj/p ~ 0. 
Normal ordering is indicated by colons. The Bessel func¬ 
tion J p of the first kind is a function of the photonic 
number operator n = a'a, where a and o' represent 
the conventional bosonic ladder operators of the res¬ 
onator. The dimensionless parameter re = Ec/Tuoq = 
h/2muj 0 is a scale related to the granularity of charges 
in the circuit via its total charging energy Eq = 2 e 2 /C, 
while Ej = Ej exp (—re/2) is a renormalized Josephson 
energy 28 . Note that experimentally Ej is easily tun¬ 
able using a superconducting quantum interference de¬ 
vice (SQUID) configuration for the JJ, from its maximum 
value down to the low percentage range. The parameter 
re is basically fixed by design and tunable in situ only 
to a limited extent, e.g., by accessing different resonator 
modes. First experiments found re ~ 0.1, while new se¬ 
tups based on highly-inductive meta-materials consisting 
of SQUID arrays 29 may already reach re ~ 0(1). 

In case that we consider the regimes of weak Josephson 
couplings going along with a low photon occupation in 
the resonator or re <C 1 such that re(n) st -C 1, the Bessel 
function can be linearized using J p (x) ~ x p /(p\2 p ). 
In lowest order, i.e., J p (-\/4rero)/n p / 2 ~ re p / 2 /p!, this 

then leads to H (p) = HAn - [(-i) p n p / 2 E*/(p\ 2)][(a t ) p + 
(—l) p a p ], which realizes special cases of systems which 
are of particular interest in the following discussion. 

For p = 1, we obtain the Hamiltonian of a driven har¬ 
monic oscillator 


= HAn + (a^ — a) , (3) 

with the driving term inducing coherent states. For p = 
2 , one has a parametric amplifier 30 

=hAn+^-((a^ 2 +a 2 }, (4) 

where the exponentiated driving term takes the form of 
a squeezing operator. For p > 2, the driving source gen¬ 
erates multi-photon processes in the resonator and thus 
nonlinear optical phenomena. 

As we will discuss in detail in the following, the full 
Hamiltonian H^ also allows to engineer V-level sys¬ 
tems by restricting the possible number of photon ex¬ 
citations in the resonator through a proper choice of the 
re parameter. For example, for p = 1 fixing re = 2, 
leads to a vanishing transition matrix element = 
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(ra I H(1) I m + 1 ) for m = l . 31 That means, the driv¬ 
ing JJ cannot cause transitions between the neighboring 
excited states 1 and 2 of the resonator and we obtain 
an effective two-level system since all inaccessible higher 
levels can be ignored (at zero temperature). 


B. Quantum master equation 


The dynamics of the reduced density operator of the 
system at zero temperature is described by a master 
equation in Lindblad form 30 


dp 

dt 


= 2 P 


h L 




+ ^ (2 apa) - np - pn ) , (5) 


which can be written in terms of the Liouvillian £. Ac¬ 
cording to the experimental realization , 6 the dissipator 
captures the effect of high-frequency modes of the electro¬ 
magnetic environment acting as a dissipative heat bath 
on the circuit which leads to photon leakage from the res¬ 
onator. The corresponding rate 7 , the inverse of the pho¬ 
ton lifetime in the resonator, is often expressed as a qual¬ 
ity factor Q = W 0 / 7 . Experimental observations 6 and 
theoretical considerations 13 show that additional low- 
frequency voltage noise affecting the JJ is weak and can 
be neglected for all observables considered here. 

In the following, we employ a decomposition, £ = £ 0 + 
3 , of the time evolution of the system 21,32,33 into two 
parts; one, describing the emission of a photon from the 
resonator in terms of a jump operator 3 defined by its 
action on the reduced density operator: 

2p = "iapa). ( 6 ) 


The remaining part £0 captures the dissipative, theoreti¬ 
cally, but deterministic system dynamics while no photon 
emission events occur. 


C. Observables for studying correlations 

Two observables will be used in the following to in¬ 
vestigate the time-independent statistics of photon emis¬ 
sion events. Their formal definition within the master- 
equation formalism is based on the Liouvillian decompo¬ 
sition £ = £ 0 + d- 

The second-order correlation function 

( 2 ), ) = (3e £T 3) st , . 

is a measure of correlations between two system jumps 
and thus two photon emission events separated by a time 
r 19,24 rjv^g no t a tion (... ) st indicates here that the sys¬ 
tem is in steady state before the first jump, i.e., (0) s t = 
Tr{Op s t} with £p st = 0. For the time-independent prob¬ 
lem here this also implies dependence on the time gap r 
only (and not individually on the two jump times). The 


time-evolution in-between the jumps is governed by the 
full Liouvillian £ allowing for additional jumps. 

The waiting-time distribution 


w(t) 


(de £ ° r 3) st 

(-j)st 


( 8 ) 


is the probability distribution for a delay t between two 
subsequent system jumps and thus two subsequent pho¬ 
ton emission events. 20,21 Now, the time-evolution is de¬ 
termined by £0 = £ — J excluding any photon emission 
events during the time delay r. 

The waiting-time distribution according to Eqs. (8) 
and (6) is sometimes referred to as a reduced form of 
a ‘full’ waiting-time distribution. The jump operator d 
is based upon ladder operators , which in the current 
context should be understood as a sum of many ‘elemen¬ 
tary transitions’, + 1 \ n ) ( n + 1|> with each tran¬ 

sition connecting only two states in Fock space. 3 is thus 
a sum over several elementary jump processes, defined as 
resulting in a certain fixed density matrix state, and so 
can not resolve the nature of the transition completely. 21 
This reflects, that from the fact that a photon leaking 
from the cavity is detected, no conclusion can be drawn 
as to which particular transition took place. 

Experimental observation of both g^ 2 \r) and w(t) is, 
of course, well established in quantum optics. In the mi¬ 
crowave regime, correlation function measurements have 
been performed by the use of linear detectors 34-36 and 
some progress has already been made for the setup con¬ 
sidered here 18 . 


III. RESULTS 

A. <r 2 (t) and w(r) at the one-photon resonance 

This first part of the results section focuses on the 
physics at the fundamental resonance, where the volt¬ 
age applied to the JJ is tuned such that the energy of 
a CP tunneling across the junction equals the energy of 
one photon at the resonator and hence each tunneling 
CP excites one photon. 


1. From the harmonic oscillator to the two-level system 

Following, we will present results for the g^(r) and 
the w(t) functions from the numerical solution of the 
Lindblad master equation (5) and discuss various ana¬ 
lytical approximations in certain regimes. Before, how¬ 
ever, we will shortly recapitulate, what to expect for the 
special cases, where the system Hamiltonian reduces to 
a driven harmonic oscillator and to a driven two-level 
system. 

For the harmonic oscillator it is well known and eas¬ 
ily shown using the Lindblad master equation, that in 
the (detuned) driven damped case its stationary solution 
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FIG. 1. (a) The second-order correlation function, g^ 2 \r), for weak driving, no = 0.2. The parameter re allows to tune the 

system from a harmonic oscillator (re = 0) to a two-level system (re = 2). (b) The second-order correlation function, 3 ®(t), 
and (c) the normalized waiting-time distribution, w(r), for the special cases of a harmonic oscillator (red lines) and a two-level 
system (blue lines). The two-level system is characterized by damped Rabi oscillations, the harmonic oscillator by a constant 
value [for g^ 2 \r)], and an exponential decay with the decay rate yno [for w(r)]. The long-time behavior of w(r) for the two-level 
system in the weak-driving regime is given by that of the harmonic oscillator, (d) The normalized waiting-time distribution, 
w(t), at re = 1.5 for various driving strengths, no- In the strong-driving limit, it shows approximately the behavior of an 
exponential decay similar to the harmonic case, however, superimposed by oscillations. 


is a coherent state |a) = exp (— |a| 2 / 2 ) Y^^Lo a<1 1 *?) /v^ 
with amplitude a fulfilling classical equations of motion 


for the two-level Hilbert space. This leads to damped 
Rabi oscillations, e.g., on resonance A = 0 , 20,37 



The corresponding density matrix reflects the defining 
property of the coherent state and is an eigenstate of the 
jump operator 0 ( 6 ) for any a. For a fulfilling the steady- 
state equation of motion (9) it is also an eigenstate of £ 
to eigenvalue 0 and thus also an eigenstate of £q. The 
time-propagation from t to t + r under £ or £ 0 then 
follows directly and yields 

Sho( t ) = !> ( 10 ) 

who(t) = 7(™)ste -7<n)atT . (11) 

Both functions clearly display the Poissonian nature of 
the photon emission events being statistically indepen¬ 
dent. In contrast to g^ 2 \r), the waiting-time distri¬ 
bution w(t) depends on a mean decay rate 7 (n) st = 
(y/nEj /h'y) 2r y/(l + 4 A 2 / 7 2 ) and thus is a function of the 
driving strength \fnEj/fry as well as the detuning A. In 
the following, we will generally use the stationary mean 
photon occupation 


"" = K (fj) (12) 

of the harmonic oscillator in case of resonance as a ref¬ 
erence measure of the driving strength even in the non¬ 
harmonic regimes. 

For a two-level system, the g^ij) and w(r) func¬ 
tions can also be straightforwardly calculated, as the time 
propagation under £ or £q can be easily explicitly solved 


3tls( t ) = 1_e " 7T cosh(/i l 7 r)+— sinh (/x l 7 r) 

wtlsO) = \n 0 e~ 1T/2 sinh 2 (^ 7 t\ 

g 2 \ z / 


J (13) 


with gi = \J (1/4 ) 2 — no and /i 2 = \J 1/4 — no- 

These results for the special cases of a harmonic oscil¬ 
lator and a two-level system, shown in the middle panel 
of Fig. 1, will in the following serve as a starting point 
to understand deviations from and the crossover between 
those special cases. The crossover is essentially governed 
by the two parameters re and rig. While the former ef¬ 
fectively changes the level structure by determining the 
relative size (and phase) of the driving matrix elements, 
the latter corresponds to the driving strength and decides 
how much of the nonlinearity of the system can actually 
be explored. 

Both the g( 2 \r) and the w(t) correlation function can 
be expressed as an expectation value of 3 with respect to 
a modified density operator p(r) = H(r) 0 p s t with H(r) 
being the respective time-evolution operator, g ^ (r) and 
w(t) can thus explicitly be calculated by means of a nu¬ 
merical implementation of the time evolution of the re¬ 
duced density operator p governed by £ and £ 0 , respec¬ 
tively. 

Numerical results are displayed in the left and right 
panels of Fig. 1, setting A = 0 for simplicity. While 
Fig. 1(a) shows how the g^ 2 \r) changes with re for a 
fixed small driving strength no = 0.2, Fig. 1(d) pictures 
the w (r) function for fixed re = 1.5 and increasing driving 
strength. 
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We easily recognize the two special cases of the har¬ 
monic oscillator and the two-level system additionally la¬ 
beled by the respective pictogram in the plot of g^(r). 
Due to the weak driving and the corresponding small 
(n) s t, only the lowest states are occupied and the system 
is close to equilibrium. g^ 2 \r) approaches here its long¬ 
time limit of 1 very soon featuring only small damped os¬ 
cillations around this value. Since the w(t) function rep¬ 
resents a normalized probability distribution, it will van¬ 
ish in the long-time limit. As in the harmonic-oscillator 
case (11), w(t) includes exponential decay terms with 
rates ~ jno- In the limit of strong driving, this exponen¬ 
tial decay, however, is superimposed by oscillations with 
frequency ~ r )y/no- For weak driving, the short-time be¬ 
havior of both g( 2 \r) and w(t) is governed by a time 
scale ~ 7 -1 . 

More features of the full numerical results will be ex¬ 
plained by the analytical solution in the weak-driving 
limit (see Sec. Ill A 2), where we will also discuss Fig. 1(a) 
further. As a preliminary conclusion, we can state, that 
our system realizes a harmonic oscillator for re = 0 and 
a two-level system for re = 2. For arbitrary values of 
re, correlations for the weakly driven system are similar 
to the harmonic case after some initial short-time devia¬ 
tions, while for stronger driving features resembling the 
two-level system appear. To understand the physical in¬ 
formation contained in the two correlation functions in 
more detail, we now turn to analytical results provid¬ 
ing exact expressions for various contributing terms and 
timescales. 


2. Approximations for weak-driving and semiclassical limits 

For the numerical results above, g ^ 2 ' 1 (t) and w(r) are 
found by calculating the time evolution over a time r of 
the stationary density matrix after a jump occurred at 
time t = 0 under the action of the Liouvillian £ [or £o 
respectively for w(t)]. For analytical results, we similarly 
consider the action of the adjoint Liouvillian £' or £q on 
system operators. Making use of the quantum regression 
theorem, 37 this yields a linear system of first-order differ¬ 
ential equations to determine <?^( T ) or w(t). However, 
in general, these systems will not close but produce a 
hierarchy of higher-order operator expressions. 

A closed system of equations will appear either if the 
number of system states is limited, or if higher-order ex¬ 
pressions drop out due to an (approximate) factorization 
into low-order terms. The number of states will thereby 
be limited if transitions to higher states are forbidden 
due to vanishing matrix elements of the driving Hamil¬ 
tonian (as discussed above in the two-level case) or small 
as in the limit of weak driving, no ~t 0. In contrast, 
higher-order expressions do not appear in the harmonic- 
oscillator case or will approximately drop out in the semi¬ 
classical regime (no 1 , re — > 0 ). 

In the regime of weak driving , no € 1, we can approxi¬ 
mately consider a three-level system and thus neglect all 



FIG. 2. The second-order correlation function, g^ 2 \r), 
in the numerical weak-driving/few-photon limit, no = 
(y/HE} / h'y) 2 —> 0, where it approaches the analytical result 
given by Eq. (14), for various values of re (with re = 0 corre¬ 
sponding to the harmonic oscillator and re = 2 to a two-level 
system). The nonlinearity for finite re becomes apparent in 
suppressed (enhanced) transitions to the second-excited state 
for re < 4 (re > 4), resulting in </ 2 ^(0) = (1 — re/2) 2 . Compare 
to Fig. 1(a) for the effects of stronger driving (see main text). 


operator contributions (al)™a ra with m > 2 or n > 2 
to find a closed system of differential equations. Solving 
this system then yields for the second-order correlation 
function the analytical expression 

<? ( 2 ) (t) = 1 + —e _7T ~ «cos (At) e -7T / 2 , (14) 


with small next-order correction terms of order no. 

Considering for the moment the resonant case, A = 0, 
the correlation takes the simple form g^ 2 \r) = (1 — 
ree -7 ' r / 2 /2) 2 so that one observes complete anti-bunching 
at times qr = 21n(re/2) for re > 2. This is indeed found 
in the numerical data of Fig. 2 for extremely weak driv¬ 
ing, while this feature tends to be absent for somewhat 
larger driving [see Fig. 1(a)]. As has already been dis¬ 
cussed elsewhere, 14 in the weak driving limit g^(r = 0) 
is directly linked to the occupation of the second excited 
state P 2 and, hence, to the probability that the resonator 
gets excited by a second CP tunneling event before it has 
relaxed to its ground state. This occupation can be found 
by a simple rate-equation description resulting in 


5 (2) (0) « 



1 

2 


Ti;. 


T, 


0,1 



(15) 


where in this regime of Coulomb-blockade a P(E )-type 
calculation 38 links the golden-rule excitation rates to the 
absolute values of the transition matrix elements of the 
driving Hamiltonian 6 . Anti-bunching, ^'(O) < 1, ob¬ 
served in the range of 0 < re < 4 is thus the effect of the 
suppression of higher-order excitations due to the nonlin¬ 
earity of the system [as reflected in the normal-ordered 
Bessel function in Eq. (2)] on the few photon level. 
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While the anti-bunching, g( 2 )( 0 ) < 1 , an indicator of 
the quantized, non-classical nature of the cavity’s electro¬ 
magnetic field can be derived from an incoherent P{E) 
description, it turns out that the dynamical features 
in (/ 2 ^(r) can not be gained such, but only in a full 
quantum-mechanical picture of coherent excitations. In 
particular, the short-time behavior g^\r) 5s g*- 2 ^(0) for 
k S; 2 can be traced back to the relative phase of certain 
transition matrix elements. 

Contrasting Fig. 2 to Fig. 1 (a) highlights the effects of 
somewhat stronger driving: 

(i) the suppression (enhancement) of g^ 2 \ 0 ) for k < 4 
(k > 4) is reduced (enlarged) as higher oscillator states 
contribute, while g^ 2 \ 0 ) = 0 for k = 2 for any driving 
strength, 

(ii) the dip occurring in g^ (r) after a short time (for 
k > 2) is reduced with g^ (r) > 0 for all times and 

(iii) oscillations with g (2 ^(r) > g > ' 2 \ oo) = 1 develop. 

Employing the corresponding approximation scheme 

for the weak-driving regime, the waiting-time distribu¬ 
tion in leading order (and again for A = 0 ) is calculated 

= e -7«or + (D e -7r _ Ke ~^/2 ( 16 ) 

77 l 0 4 

with small next-order correction terms of order no ■ Com¬ 
paring to Eq. ( 14 ), the close relation between g^(r) and 
w(t) for weak driving is apparent. 

Apart from their different normalization, we observe 
that the second and third terms of both expressions 
with the decay rates 7 and 7/2 [and thus the short- 
time behavior of g^{r) and w(t)] coincide. This hap¬ 
pens since the underlying systems of differential equa¬ 
tions based on £ and £ 0 differ in higher-order operator 
expressions which do not contribute in leading order in 
no to these two terms (but show solely an effect on the 
first term). Nonetheless, for long times w{r) finally shows 
a slow exponential decay on a time scale 1/(7770), while 
g^(r) —► 1 in this limit. 

Notably, the time scale (7/2) -1 which impacts the 
time-dependent behavior of both g^ (r) and w(t) in the 
strongly overdamped case can be traced back to the de¬ 
cay rate of coherences, i.e., off-diagonal entries of the 
density matrix. This, in turn, means, that in contrast to 
g( 2 \r = 0 ) ( 15 ) the time dependence of the correlation 
functions can not be properly described within a P(E)-like 
rate-equation approach even in the weak-driving limit. 
This also highlights the profound difference of the pho¬ 
tonic two-level case in comparison to the case of correla¬ 
tion functions of (spinless) electrons transported through 
a single site, where coherences do not matter 21,24 . 

In the semiclassical regime , the operators can be 
replaced by a/*) + Sa"\ with Sa^ describing quantum 
fluctuations around the classical solution of the complex 
oscillation amplitude a/*) = (a^). In the semiclassical 
regime, fluctuations will then be small compared to the 
mean amplitude and operator contributions beyond lin¬ 
ear order can be neglected. 
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FIG. 3. The second-order correlation function, g^ 2 \r), 
for strong driving, no = 18 2 . At special values of k (dot¬ 
ted lines) certain transition matrix elements vanish, reduc¬ 
ing the number of involved cavity states. This yields pro¬ 
nounced (damped) oscillations in </ 2) (r) in a certain k inter¬ 
val around these values. Lines below k, = 2 can be traced 
back to the first zero at 3.83 of the Bessel function Ji, e.g., 
T2,3 oc (2 j : Ji(y/ 4 m)jffn : | 2) = 0 for k = 1.27, the line at 
k = 3.31 to the second zero. 


Such a semiclassical approach proved to be very fruitful 
in calculating stationary expectation values revealing, for 
instance, bifurcations between different solutions. 12,16 

In principle, one can proceed similarly here and calcu¬ 
late c/( 2 )(t) and w{r) in the semiclassical limit. It turns 
out, however, that in the regime, where the semiclassi¬ 
cal approximation is valid, the resulting behavior of both 
functions is hardly distinguishable from that of a har¬ 
monic oscillator: g^ is basically 1 [Eq. (10)] with tiny 
corrections of order l/|a| 2 and w(r) is dominated by a 
strong exponential decay [Eq. (11)] due to a large (n) st , 
which completely obscures small (quantum) corrections 
occurring in the long-time limit. 


3. Effective N-level systems 

Above, we explained, how for the special value of k = 2 
the system is reduced to a two-level system. Analyz¬ 
ing the roots of higher transition matrix elements be¬ 
tween neighboring states reveals the occurrence of other 
few-level systems for specific choices of n. (See Ref. 39, 
where AT- level systems are realized by dynamically block¬ 
ing transitions with an additional external drive.) For 
example, the conditions X 2 i3 = 0 valid for k = 3 ± \/3 
and T 3j4 = 0 valid for n ss 0.94, 3.31, or 7.76 lead to 
effective three- and four-level systems, respectively. 

Note that the vanishing of transition matrix element 
Tm, m + 1 oc (777 | : J i(V ’inn) / y/n : | to) corresponds to ze¬ 
ros of the Bessel function J\ for large to. Due to normal 
ordering, however, for small to th e roots are shifted, e.g., 
for to = 1 we have (l | J\(\j 4ren )/yjn \ l) = Ji(V4k) = 
0 for k = 3.67 without the normal ordering instead of 
k = 2 . 
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FIG. 4. The second-order correlation function, g (2 \r), as a 
function of detuning, A, at re = 1 for weak driving, no = 0.2. 
Oscillations with frequency / ss A/2n are observable with an 
amplitude, which increases with increasing detuning. Com¬ 
pare to the analytical result given by Eq. (14) for the effects 
occurring in the weak-driving limit, no = {y/nEj /h'y) 2 —> 0. 


Figure 3 displays the second-order correlation function 
as a contour plot in the regime of very strong driving 
(no = 18 2 ) for re £ [0, 4], We immediately recognize the 
prominent Rabi oscillations in the case of a two-level sys¬ 
tem, see Fig. 1(b). However, also TV-level systems feature 
Rabi-like oscillations indicating those special values of re 
where higher transitions are suppressed. The amplitude 
as well as the frequency of these oscillations decrease with 
an increasing number of remaining levels. Notably, the 
Josephson-driven resonator setup thus allows by properly 
adjusting the re parameter to turn the harmonic oscilla¬ 
tor into a highly nonlinear TV-level system. If re could be 
varied in situ (for example by an external magnetic flux 
in case of a SQUID array), one could even switch from 
TV to TV' level systems on short-time scales. 


4- Detuning from the one-photon resonance 

We drop now the restriction of A = 0 to study the 
effect of detuning on the system. Considering the special 
case of a harmonic oscillator, the analytical results for 
both g( 2 )(r) [Eq. (10)] and w(t) [Eq. (11) are still valid 
without limitations but it should be noted that (n ) s t is 
modified by a factor of 1/[1 +4(A/7) 2 ], thus leading to a 
much slower decay for high detuning. Leaving this special 
case and going to higher re, the g^ and w{r) functions 
are both superimposed by oscillations. In the regime of 
weak driving, the behavior of g^ (r) is directly given 
by Eq. (14) revealing oscillations with frequency w = A 
and high amplitudes for strong detuning, whereas g^(0) 
remains unchanged. The numerical results of g^ (r) for 
somewhat stronger driving (no = 0.2) in Fig. 4 show 
small differences to the analytical results, only. Note 
that the far off-resonant two-level system also shows Rabi 
oscillations with frequency A for strong driving. 


B. Relations between g^ 2 \r) and w(r) 

In order to investigate the statistics of the photon emis¬ 
sion events in various regimes of our system, we have so 
far discretionarily elected to discuss at times the second- 
order correlation function g^ (r) or the waiting-time dis¬ 
tribution w(t). Obviously, this immediately throws up 
the question, to what extent the information provided by 
gl 2 -*(r) and w(t) differ or complement each other. 

Apparently, both tools oftentimes provide very simi¬ 
lar information and are closely linked. Consider, for in¬ 
stance, that their definition [Eqs. (7) and (8)] directly 
yields the equal-time relation w(0) = 7 (n) st g (2 ^( 0 ), and 
that we have observed similar Rabi-type oscillations in 
both functions, when certain transitions between excited 
neighboring states are suppressed. 

In general, however, there is no one-to-one correspon¬ 
dence between g( 2 )(r) and w(t). In the limit, where each 
relaxation process resets the system to a known state, 
a situation described in this context by renewal theory, 
such an equivalence exists though. To see that, one works 
with characteristic second-order Fano factors (see follow¬ 
ing) for the respective quantities. 


1. Renewal theory 


A Poissonian process, where events are occurring inde¬ 
pendently, can alternatively be considered as character¬ 
ized by the fact that after each event the probability for 
the subsequent event is given by an identical, indepen¬ 
dent exponential probability distribution. Renewal the¬ 
ory generalizes this to so-called renewal processes charac¬ 
terized by an identical, independent, but otherwise arbi¬ 
trary probability distribution. 40 Hence, each event resets 
the system to the very same state and thus renews it. 

For such a renewal process, the two functions g^ 2 \r) 
and w(t) can be directly related (in Laplace space), 24 and 
thus provide identical information. From that relation a 
somewhat simpler comparison restricted to low-order mo¬ 
ments follows. 22,23,25 Indicating the relative strength of 
second-order fluctuations, a single, dimensionless num¬ 
ber, the Fano factor is used. 

The standard Fano factor, F fcs , characterizing the 
statistics of the photon flux J p h leaking from the cavity, 
can in the spirit of full counting statistics (FCS) be de¬ 
fined from the statistic of the number of photons leaked 
during an accumulation time, Nt = dr / p h(r) as 


pFCS 


(N£) - (Nt) 2 
(N t ) 


T—>oo 


(17) 


Following Ref. 19 it can also be gained from g^ 2 \r), 


pFCS 


l + 2(/ ph W dr g^ 2 \r) — 1 


S(u = 0) 

2</ph> 


(18) 














FIG. 5. In case that renewal theory applies to the system, g (2 \r) and w(r) provide identical information and the two Fano 
factors deduced from them are identical: _F fcs = F WTD . The measure |j? wtd /F fcs — 1| (a) highlights the deviations of the 
system from renewal theory, which vanish at k = 0 (harmonic oscillator), k = 2 (two-level system), and no < 1 (weak-driving 
regime, where nonlinearities are not crucial). The dotted lines refer to two cross sections (b) at weak and strong driving, 
revealing sub-Poissonian distributions of emitted photons. 


where the second line highlights an alternative access to 
this quantity via the zero-frequency noise, S(u> = 0), of 
the photon current. 

An alternative Fano-type factor, F WTD , can be defined 
from the first two cumulants of the waiting-time distri¬ 
bution, 


Then, for any renewal process F fgs = F wtf> holds with 
F fgs = 1 = F wtd for the Poissonian case. 41 For our JJ- 
resonator system it is obvious that in the two-level case, 
k = 2, the system is reset to the very same state, namely 
|0)(0|, by any jump process and, hence, the photon emis¬ 
sion process constitutes a renewal process. 

For the harmonic oscillator, realized for k —> 0 or in 
the weak-driving limit, the jump operator acting on an 
arbitrary system density matrix does not actually reset 
the system to a single state. It does so, however, for 
any density matrix appearing during the time evolutions 
within the definitions [Eqs. (7) and (8)] of g^ 2 \r) and 
w(t) , since the time evolutions start from the stationary 
state of the harmonic oscillator, an eigenstate of 3 and 
of both £ and £o- Thus, for our purposes the harmonic 
oscillator behaves as a single-reset system and undergoes 
a renewal process. 


difference vanishes, which can be observed in Fig. 5(a) 
in a small region around the special cases of n = 0 and 
k = 2 for arbitrary driving strength. It also holds true 
for arbitrary n in the regime of weak driving, where the 
dynamics of the system are hardly influenced by the non- 
linearities of the Bessel function. Strong discrepancies 
from a renewal character appear in the classical regime 
at small but finite k and strong driving no- The shape 
of that feature, visible in Fig. 5(a), reflects that nonlin¬ 
earities appear, once t he a rgument of the Bessel function 
becomes of order 1, \/4 nn ~ 1. 


The two cross sections at no = 0.5 and 4.0 in Fig. 5(b) 
showcase two facts. First, we observe, that there may be 
additional parameter values with F fgs = F WTD while 
renewal theory is not expected to be valid. Indeed, for 
these values, one can check, that the full time depen¬ 
dence of g^ 2 \r) and w(t) does not give identical infor¬ 
mation, but differences only become apparent in higher 
moments. Secondly, note that the sub-Poissonian dis¬ 
tribution of emitted photons in the long-time limit with 
F fcs < 1 for k > 4 does not correspond to anti-bunched 
photons, since, in fact, g^( 0) > g^ 2 \r) at all times for 
these values of k (and weak to moderate driving) [see 
Figs. 1(a) and 2], 


2. JJ-resonator system and its renewal character 

In order to capture the system’s deviations from re¬ 
newal character, the two different definitions of Fano 
factors, F fcs and _F WTD , are compared in Fig. 5. If 
the two Fano factors differ, g^ 2 \r) and w(t) are not di¬ 
rectly linked and provide, in principle, different informa¬ 
tion about the system. In case of a renewal process, the 


Furthermore, we observe F fgs > 1 for strong driving, 
e.g., choosing n 0 = 18 2 , at n = 3 — \/3 corresponding to 
an effective three-level system, although we clearly expect 
here an anti-bunching behavior since g^(0) < g^ 2 \r) 
at all times (see also Fig. 3). These findings are in line 
with recent results 24,42-44 emphasizing that Eq. (18) does 
not in general imply a one-to-one correspondence of sub- 
/super-Poissonian FCS and short-time (anti-)/bunching. 
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FIG. 6. The waiting-time distribution (solid line), w(t), for 
the two-photon resonance in the special case of a parametric 
amplifier, k = 0.001, nQ Ph = 0.02 <g 1, which approximately 
coincides with the analytical result given by Eq. (21). The 
short- and long-time decays are indicated by the analytical 
results (dashed lines), corresponding to the probability of de¬ 
tecting the first or second photon of an emitted pair. 


C. g and w(r) at the two-photon resonance 

So far, only small detuning around the resonance con¬ 
dition uij = cuq for one-photon processes have been inves¬ 
tigated. Tuning the bias voltage, other than this funda¬ 
mental resonance can be accessed. In particular, we now 
turn to the simplest higher-order creation processes, the 
two-photon process with resonance condition 2ui 0 ss wj. 
We start our considerations on the basis of the Hamil¬ 
tonian given in Eq. (2) with p = 2, which is in a frame 
rotating with wj/2 and within a rotating-wave approxi¬ 
mation. For most of the following, we consider the regime 
of weak driving and consequently small occupation (n) s t, 
where the Hamiltonian simplifies to that of the paramet¬ 
ric amplifier, in Eq. (4). 

In this weak-driving regime, incoherent CP transport 
across the JJ can be observed, i.e., the events of CP tun¬ 
neling processes are statistically independent and follow 
a Poissonian probability distribution. This is based on 
the fact that photon relaxation at the resonator occurs 
sufficiently fast compared to CP tunneling. Therefore, 
we always observe a pair of photons followed by a long 
waiting time until the next emission event occurs. 

By elementary considerations following this picture, 
the two different Fano factors _F fcs ^ F WTD are found 
to differ, since renewal theory does not hold. The ob¬ 
servation that for well-separated bunches of p parti¬ 
cles the Fano factor gives the bunch size, F fcs = p , 
is well known (see, e.g., Refs. 45-48). Formally, one 
may argue that in the long-time limit, the number of 
photons Nt = / q T drI p h(r) is well approximated by a 
coarse-graining over typical time scales of the bunch 
by Iph (t) « plcp(r), where the CP (number) current 
Jcp captures the underlying (Poissonian) statistics of the 
bunches. The Fano factor then immediately follows from 


Eq. (17), F fcs = ^-Fq p S = p with corrections deter¬ 
mined by the separation of the typical time-scales of the 
bunch length and the time between bunches. 

The first and second moments of the waiting-time dis¬ 
tribution, and hence F' WTD , are found averaging over 
the equally probable cases, that the first or the second 
photon of a bunch (for p = 2) is observed at t = 0. 
This gives (r ph ) = [ 0 ( 7 _1 ) + (r C p)]/2 and (r^ h ) = 
Ph- 1 ) + ( t Cp)]/ 2 = 4(T ph ) 2 with corrections on the 
order of the bunch length, 7 -1 , resulting in F WTD = 3 
(for general p , one similarly finds F' WTD = 2 p — 1 ). 

Going beyond these preliminary considerations, ana¬ 
lytical expressions for the full time-dependence of g ^ (r) 
and w(t) can again be found in the weak-driving limit 
(restricted to the resonant case, ujq = wj/2, again). 

The instantaneous limit, 

3 (2) (0) = -4ph+2 (20) 

zn 0 


with next correction terms of the order of the stationary 
mean occupation number ng ph = (kEj /fi. 7) 2 /2 (in lead¬ 
ing order), which again measures the driving strength, 
has been extensively discussed elsewhere. 11 ’ 13 ’ 14 "We want 
to emphasize, that the l/2ng ph -divergence of photonic 
correlations in the weak-driving limit due to the two- 
photon creation process can again be understood on the 
basis of excitation and decay rates for double and single 
occupancy of the cavity, while the time-dependence of 
<? (2) (r) can not be gained correctly in this simple man¬ 
ner. Note also, that depending on the quality factor 
of the cavity, two-photon processes may dominate the 
photonic correlations even around the fundamental res¬ 
onance, where they compete with the resonant single¬ 
photon processes. 14 

The waiting-time distribution, shown in Fig. 6, takes 
the form 


w(t) 

7 


2ph 


3 -7T . 


^0 P V/2 


13n, 


2ph 


27T 


( 21 ) 


with both coefficients and decay rates in leading order 
in nl ph . The expectation value (r) for the waiting time 
is dominated by the first and the second contributions, 
while the third one does hardly contribute for weak driv¬ 
ing when n 2ph -C 1. Physically, the first contribution, 
describing the case that the first photon of a bunch is 
observed at r = 0, determines the waiting-time for the 
second photon within the same bunch. Instead, the sec¬ 
ond contribution is related to the waiting-time for the 
next bunch to arrive, when the second photon was ob¬ 
served initially. 

Using the analytical results of </ 2 H r ) and w(t) in 
conjunction with Eqs. (18) and (19) results in F fcs = 
2 + 7ng ph /2 and F WTD = 3 + 3nQ Ph confirming our above 
considerations. 








10 


IV. CONCLUSIONS AND OUTLOOK 

We have analyzed the time-resolved statistics of pho¬ 
ton emission events of an electromagnetic oscillator cou¬ 
pled to a voltage-biased Josephson junction by means 
of the second-order correlation function </ 2 )(r) and the 
waiting-time distribution w(t). Starting from a time- 
independent Hamiltonian in rotating-wave approxima¬ 
tion and the quantum master equation in Lindblad form, 
we have obtained numerical results as well as analytical 
expressions on the basis of perturbative approaches in 
the weak-driving and the semiclassical regime. The pho¬ 
tonic statistics are basically governed by the parameters 
Ej/fvy measuring the driving strength and n describing 
the impact of charge quantization and determining mag¬ 
nitude and phase of the transition matrix elements of the 
drive. 

Changing parameters allows the realization of various 
simple systems and thus the observation of the pho¬ 
ton statistics of a driven harmonic oscillator, a two- 
level system or a parametric amplifier, all within the 
voltage-biased JJ-resonator compound. Increasing the 
impedance and hence k. either by design or in situ by 
use of highly inductive metamaterials, the crossover from 
the harmonic oscillator (k —> 0) to the two-level system 
(k = 2) can be studied. For other special values of k , 
further iV-level systems can be implemented, exhibiting 
for instance Rabi-type oscillations in g^ 2 \r) and w[r) for 
strong driving. 

While the scenarios above occur at the fundamental 
resonance where each tunneling Cooper pair creates one 
photon, higher resonances can be accessed when the volt¬ 
age bias provides the energy necessary to create several 
photons. Then, for example, at the two-photon reso¬ 
nance, parametric-amplifier physics can be investigated 
in the n —» 0 limit. 

Photon statistics will show highly characteristic sig¬ 
natures in these different scenarios: most dramatically, 


single-photon creation and complete anti-bunching are 
observed for the two-level system, photon-pair creation, 
and bunching for the parametric amplifier. To further 
study the full time-resolved photon statistics, we ex¬ 
plored the second-order correlation function g^ 2 \r), the 
waiting-time distribution w(t), and Fano-type factors ex¬ 
tracted from them. The rich variety of strongly corre¬ 
lated and non-classical states of light, observable even 
for weak driving, all stem in essence from the inher¬ 
ent nonlinearity of the Josephson junction appearing as 
a normal-ordered Bessel function in the effective RWA 
Hamiltonian. 

Extending our study to several modes within one or 
several cavities will be of particular interest for further 
research. If the energy of a CP transfer equals the sum of 
the energies of two photons, entangled photon pairs can 
be created. The statistics of these photon pairs have so 
far only been investigated by means of the second-order 
correlation function and Fano factor, 15,16 however, not 
by means of the waiting-time distribution. 

Furthermore, a system consisting of JJ and resonator 
makes it generally possible to observe both photonic 
and current statistics (Josephson photonics). Of course, 
the dynamics of the charge transfer processes can 
indirectly be investigated by observing the created 
light. However, it would be natural to also directly 
consider the waiting-time distribution for subsequent 
charge transfer processes at the JJ or even a mixed 
waiting-time distribution relating CP transport statistics 
to photonic emission events. 
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